% inverse image and mu given as an array of Numcolors X RGB
function [strip all] = wcd2(im, mu)

% Take the Red component
[height width channels] = size(im);

sim = imresize(im,0.2);

numofcolors = size(mu,1);

V = [];
C = [];
for colorIdx =1:3
    v = [];    
    c = sim(:,:,colorIdx);
    
    npix = size(c,1) * size(c,2);
    
    for cidx = 1:numofcolors
        vi = speye(npix);
        vi = vi * mu(cidx,colorIdx);

        v = [v,  vi]; 
    end
    
    V = [ V ; v ];
    C = [ C ; c(:) ];
    
end


%t = reshape(C, size(im,1), size(im,2), 3);
%imshow(t);

x = lsqnonneg(V, double(C));
%x = V \ double(C);

tnr = V * x;
tnr = reshape(tnr, size(sim,1), size(sim,2), 3);

rim = imresize(tnr, [height width] );
t = [im rim];
figure(1);
all = [t; 1-t];
imshow(all);


s = [];
strip = [];
for cidx = 1:numofcolors
    delta = (cidx - 1) * npix;
    selector = delta+1:delta+npix;
    nr = V(:, selector) * x(selector);
    nr = reshape(nr, size(sim,1), size(sim,2),3);
    
    %hist(nr);
    figure(2);
    t = 1 - nr;
    t = imresize(t, [height width] );
    
    c = mu(cidx,:);
    c = reshape(c, [1 1 3]);
    c = imresize(c, [height width] );
    
    t = [t ; 1-c];
    
    if(isempty(s))
        s = nr;
    else
        s = s + nr;
    end
    
    strip = [strip t];
end


imshow(1 - s);





